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Abstract 

The J-PET scanner, which allows for single bed imaging of the whole human 
body, is currently under development at the Jagiellonian University. The dis¬ 
cussed detector offers improvement of the Time of Flight (TOE) resolution due 
to the use of fast plastic scintillators and dedicated electronics allowing for sam¬ 
pling in the voltage domain of signals with durations of few nanoseconds. In this 
paper we show that recovery of the whole signal, based on only a few samples, 
is possible. In order to do that, we incorporate the training signals into the 
Tikhonov regularization framework and we perform the Principal Component 
Analysis decomposition, which is well known for its compaction properties. The 
method yields a simple closed form analytical solution that does not require iter¬ 
ative processing. Moreover, from the Bayes theory the properties of regularized 
solution, especially its covariance matrix, may be easily derived. This is the key 
to introduce and prove the formula for calculations of the signal recovery error. 
In this paper we show that an average recovery error is approximately inversely 
proportional to the number of acquired samples. 
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1. Introduction 

Positron Emission Tomography (PET) [TJ[^[3] is currently one of the most 
prominent and promising techniques in the field of medical imaging. It plays a 
unique role both in medical diagnostics and in monitoring effects of therapy, in 
particular in oncology, cardiology and neurology. Therefore, notable efforts are 
devoted to improve this imaging technique. The best way so far is to determine 
the annihilation point along the Line of Response (LOR) based on measurement 
of the time difference between the arrival of the gamma quanta at the detectors, 
referred to as the Time of Elight (TOE) difference [H |S]. As it was shown in 
Ref. [6], even with the TOE resolution of about 400 ps that is achievable with 
non-organic crystals, a signal-to-noise ratio can be improved substantially in 
reconstruction of clinical PET images. 

In the articles laiaiaiin], a new concept of the TOE-PET scanner, named 
J-PET, was introduced. The J-PET detector offers improvement of the TOE 
resolution due to the use of fast plastic scintillators. A single detection unit of 
the newly proposed TOE-PET detector is built out of a long scintillator strip. 
Light pulses produced in the strip propagate to its edges where they are con¬ 
verted via photomultipliers into electric signals. There are two main reasons 
why the TOE resolution may be improved in J-PET scanner: i) a very short 
rise-time and duration of the signals and ii) a relation between the shape and 
amplitude of the signals and the hit position. The latter feature usually dis¬ 
torts the time resolution but, when the waveform of the signal is registered, 
the information about a change of the shape with the position may increase 
the position resolution and indirectly improve also the resolution of the time 
determination m- However, to probe the signals, with duration times of few 
nanoseconds, a sampling time of order of picoseconds is required. This can be 
done well with the oscilloscopes during the laboratory studies on the proto¬ 
type, but in the final multimodular devices with hundreds of photomultipliers, 
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probing with oscilloscopes is not feasible [I1II31. Therefore, sampling in the 
voltage domain using a predefined number of voltage levels is needed. An elec¬ 
tronic system for probing these signals in a voltage domain was developed and 
successfully tested [la¬ 
in recent papers we have investigated the performance of a single 

unit of a J-PET scanner. Sampling in the voltage domain at four thresholds 
was simulated and each pair of waveforms was represented by a 15-dimensional 
vector holding information about the relative time values of a signal’s arrival 
at both scintillator ends [T^. In that scenario, the spatial and time resolutions 
of the hit position and event time for annihilation quanta were determined to 
be 1.05 cm and 80 ps (cr), respectively. It is evident that the spatial and time 
resolutions can be further improved primarily by an increase in the number of 
threshold levels, as was also concluded e.g. in article m- However, the number 
of channels in the electronic devices is a very important factor in determining 
the cost of the PET scanner. Therefore, the question arises: is it possible to 
recover the whole signal based on only a few samples? Equivalently, how many 
threshold levels have to be applied to achieve a reasonable estimation error? 

In this article we propose a novel signal recovery scheme based on ideas from 
the Tikhonov regularization mill] and Compressive Sensing methods 

that is compatible with the signal processing scenario in J-PET devices. We 
investigate the quality of signal recovery based on the scheme with a single 
scintillator strip module introduced in Ref. mm- The two most important 
aspects of our work involve i) a development of fast recovery algorithms and 
ii) a statistical analysis of an error level. In practice the algorithm needs to 
work in real-time scenarios: during a single PET examination more than 10 
million signals are acquired in just 10-15 minutes. Moreover, only results for 
realistic scenarios with noisy measurements are considered. In particular, as 
was mentioned, the most important part of our investigations is to determine a 
dependence of the signal recovery error on the number of samples taken in the 
voltage domain. In this paper the formula for calculations of the recovery error 
will be introduced and proven. 
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This article is organized as follows. We will define the problem of signal 
recovery and show briefly the Tikhonov regularization and compressive sensing 
methods in Sec. 2. In the last part of this section we will introduce the theorem 
enabling the determination of the signal recovery error as a function of the 
number of samples. The experimental setup of the simplified PET device with 
a single scintillator strip that enables us to acquire real signals as well as the 
results of their analysis are presented in Sec. 3. A detailed analysis of the 
experimental characteristic of signal recovery error as a function of the number 
of samples, as well as the explanation of the specificity of the signal recovery 
method in the application to the J-PET measurement are provided in Sec. 3.2. 
In Sec. 3.3 we have discussed the limitations of the method of signal recovery. 
In particular, we have presented how the quality of the information needed to 
recover the signals, and therefore to estimate the recovery error, vary with size 
of training set of fully acquired signals. Moreover, we have demonstrated that 
using the recovered waveform of the signals, instead of samples at four voltage 
levels alone, improves the spatial resolution of the hit position reconstruction. 
A detailed description of this study is given in Sec. 3.4. The conclusions and 
directions for future work are presented in Sec. 4. 

2. Materials and methods 

2.1. Problem definition 

We wish to recover a finite signal G TZ^ in a situation where the number 
M of available samples, denoted as measurement y G TZ^, is much smaller 
than the signal dimension N (y is sampled on some partial subset fl, where 
the cardinality |f2| = M). In the compressive sensing (CS) method [501 [H], a 
sparse expansion G TZ^ of signal y^, evaluated via linear and orthonormal 
transformation y^ = Ax^, is considered. In the following we assume we are given 
a contaminated measurement y and then one may write: y = Aq^x'^ + e; where 
An is a M X N matrix modeling the sampling system, constructed from M rows 
of matrix A that corresponds to the indexes of y described in the subset fl. 
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and e is an error term. Therefore, during the recovery process the information 
about the measurement y may be included in the form of the linear system of 
equations: 

y = Aqx. (1) 

It should be stressed that in the case of presence of noise, represented by signal 
e, instead of an exact recovery of signal we will consider the solution x and 
by the analogy, instead of signal y'^ we will consider the solution y. Before we 
look in more detail we may state that the evaluation of y requires two steps: i) 
recovery of the sparse expansion x and ii) calculation of y based on the x. The 
first step of the procedure is crucial. 

In the situation where an exact solution cannot be found, CS method provide 
an attempt to recover x by solving optimization problem of the form 

a; = argmin ||a:||i such that ||j/— ^na ;||2 < e 

where e is the size of the error term e. The li minimization approach provides 
a powerful framework for recovering sparse signals. Moreover, the use of li 
minimization leads to a convex optimization problems for which there exist a 
variety of greedy approaches like Orthogonal Matching Pursuit [52] or Basis 
Pursuit |23|. Other insights provided by CS are related to the construction of 
measurement matrices (Aq) that satisfy the Restricted Isometry Property [2U 
[25] . For an extensive review of CS the reader is referred to Refs. [SolETlIMlES]. 

We will incorporate from the CS framework to our scheme only the idea of 
conducting the experiment, formulated by a linear system of equations as given 
explicitly in Eq. Q. The problem formulated by Eq. Q alone is essentially 
underdetermined, and is so-called ill-posed [26]. As in the case of the CS method, 
it is necessary to incorporate further assumptions or information about the 
desired solution in order to stabilize the problem. As an alternative to the CS 
theory, one may use the regularization methods nilmiETlEH] . The Tikhonov 
regularization (TR) method [THl[in] is the most suitable for our problem. Here, 
the idea is to define the regularized solution x as the minimizer of the following 
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expression: 

X = argmin{(?/ - A^xY"R~^{y - Ana;) + (a: - YfP~^{x - y.)}- (2) 

In Eq. Q, both signals y and x are assumed to be given with multivariate 
normal (MVN) distributions, 

y ^ N{Anx,R), (3) 

a: ~ A/'(/i,P), (4) 

where Ana; and R are the mean value and covariance matrix of a measured 
signal y, respectively, y and P are the mean value and covariance matrix of 
a prior distribution of x^, respectively. The covariance matrix R in Eq. (|^ 
is diagonal with the values on the diagonal equal to the measurement error 
variances (as explained in Sec. 2.4). With the introduction of the second term 
to the optimization problem in Eq. ([^ , an additional information from a training 
set of fully sampled signals y^ is provided. The prior distribution of sparse 
representation a;° (see Eq. @) is evaluated based on the linear transformation 
of the training set of signals y^ by using the Principal Component Analysis 
(PCA) decomposition [53]. Thus, in order to find the sparse representation x of 
a given measurement y, as a solution of Eq. ([^, one needs to specify hrst the 
prior distribution of x^. 

Beside the advantage of including the additional information from training 
signals, a further benefit of the TR approach is that the problem in Eq. (2) 
has an optimal solution which can be determined explicitly. In Sec. 2.2 we 
will evaluate the orthonormal matrix A, as well as the parameters of the prior 
distribution of signal a;°, see Eq. Q, via the PCA decomposition of training 
signals j/°. It should be stressed that these parameters are calculated only once, 
at the preparation stage of the procedure. Thus, the same matrix A, P and 
vector /r, are used to recover a signal x for each measurement y. An example 
of the idea of using the PCA decomposition of a training data set in a similar 
problem may be found in Ref. m- In Sec. 2.3 the solution of the TR formula 
described in Eq. Q, as well as its properties, will be provided. Finally in Sec. 2.4 
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we will introduce the theorem enabling the determination of the signal recovery 
error as a function of the number of samples. 

2.2. Principal Component Analysis 

PCA is a statistical study, based on the orthogonal transformation, to con¬ 
vert a set of signals into a set of linearly independent variables, such that the 
variance of the projected data is maximized. For the training data matrix of 
fully sampled signals 

^ = [(2/°!) - m)\iyl2) - HI-Ky^L) - w)], (5) 

where m is the mean of the aligned L training signals y^, the PCA coordinates 
A = [a;^^)|a;^ 2 )l-| 2 ;(L)] are given by: 

A = A^Y. (6) 

The matrix A = [a(i)|a( 2 )|...|a(Ar)] in Eq. ^ is calculated in such a way that the 
projection of the data matrix Y with successive basis vectors a(i), a( 2 ),..., a(Ar) 
inherits the greatest possible variance in the data set Y. Thus, the first basis 
vector has to satisfy: 

«(!) = argmax||a'^r||2, 

where ||a ||2 = 1 (the orthonormality is restricted). The component can be 
found by subtracting the first k — 1 principal components from data set Y: 

k-l 

Efe = Y 

1=1 

and then finding the basis vector which extracts the maximum variance from 
this new data matrix 

a(k) = argmax||a^Tfc||2, 

where ||a ||2 = 1. 

In the case discussed in this paper, the matrix A is evaluated based on the 
PCA decomposition of the training set of signals and therefore the parameters 
of the MVN distribution of {fi, P) are estimated based on data matrix A, 
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constructed according to Eq. Q. The empirical covariance matrix P of data 
set X may be evaluated as: 


P = E[X ■ X^]. 


(7) 


The covariance matrix P is diagonal, with values sorted in non-increasing order. 


Since the mean of the signals in data set Y is equal to 0, see Eq. ([^, the mean 


fi = 0. 


2.3. Tikhonov regularization 

In the previous section we have shown how the prior information from a 
training set of signals i.e. the orthonormal matrix A, and the parameters 
of the prior distribution of signal may be introduced to the TR framework. 


In this section we will derive a sparse solution x of Eq. ([^ , and its covariance 


matrix, denoted hereafter as S, for a particular measurement y, based on the 
TR assumptions [HI [19]. The posterior probability density function (pdf) of the 
signal X conditional on measurement y, namely p{x\y), can be computed after 
combining the prior distribution of x, p{x), likelihood of measurement p{y\x), 
and p{y) via the well-known Bayesian rule: 


r I ^ p{x)-p{y\x) 

P(x\y) = . ' 


( 8 ) 


To describe the MVN distribution in Eq. ([^ we will use the following notion: 



where z is an iV—dimensional variable with mean value u and covariance matrix 


Q. Hence, the marginal and conditional densities of x and y from Eq. (§ are 
given as follows: 


p{x) =N'{x\p,P), 
piy\x) =J\f{y\Anx,R) 
p{y) = a, 


(9) 

( 10 ) 

( 11 ) 

( 12 ) 


p{x\y) =M{x\x,S). 
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Equations Q and (10) result directly from the previously described Eq. (|^ 
and ([^, respectively. Equation © shows that the probability p{y) is inde¬ 
pendent of x, and therefore serves as a normalization constant. The posterior 


probability in Eq. (12) can be described exclusively by its first two moments 


{x, S) because a Gaussian pdf is self-conjugate and the pdfs on the right hand 
side of Eq. (|^ are Gaussian. After some simple calculations the equations for 
X and S are given by m- 


X — {P ^p + AqR ^y) ■ {P ^ + AqR 
S = {p-^ + A^R-^An)-\ 


(13) 

(14) 


It is worth noting that the solutions in Eq. (13) and (14) are analogous to 
Kalman filter update equations (cf. Refs. [311[33]). It can be easily shown that 
X is not only the minimum mean square error (MSE) estimator (see Eq. (§) 
but also the maximum a posteriori (MAP) estimator, that is 

X = argmax{p(x|j/)}. 

It should be stressed that all the information from the training set of signals 
(matrix A, P and vector y) and from the oscilloscope specification (matrix 
R) are evaluated only once, at the preparation stage. Thus, the sparse signal x 


may be found, according to Eq. (13), as a linear combination of the previously 


defined parameters and a given measurement y. However, the evaluation of the 


covariance matrix S, according to Eq. (14), does not require the information 


about the measurement y, and may be provided at the preparation stage. This 
fact opens a possibility for an estimation of the theoretical value of the recovery 
error. This idea will be presented in the next section. 

2.4- Analysis of the signal recovery error 

As mentioned at the beginning of the Sec. 2, the evaluation of the recovered 
signal y requires two steps: i) recovery of the compact representation x via 


Eq. (13) and ii) calculation of y as the solution 

ij = Ax + m, 


(15) 
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where m and A are derived by PCA decomposition. One of the benefit of 
using the TR approach is that it provides an easy way to obtain the error 
term of the recovered signal y. We assume for the sake of simplicity that m 


in Eq. (15) is known exactly. Since the matrix A is orthonormal, we have 
||y~ 2/°||2 = 11^^ — a ;°||27 and therefore we may focus on the recovered signal x 


error. 

In multivariate statistics, the trace of the covariance matrix is considered as 
the total variance. We will denote the trace of covariance matrix S' as cr^. It 
is worth noting that cr^ is the mean value of the recovery error squared norm 
||x — Let P{k) be the diagonal element of covariance matrix P (see 

Eq. 0). Find the smallest value D, and largest value r (with constraints D > 0 
and T > 0) such that for each 1 < k < N: 


P{k) < D ■ 


(16) 


From Eq. (16) one may see that t controls the decrease rate of P(k)\ the greater 
r, the faster the decreasing of P(k) and better the compressibility of signal x. 
The characteristics D and r of the prior distribution of signal x and a standard 
deviation of noise (cr) enable us to provide the formula for average value of the 
recovery error tr^. For this purpose we formulate the following theorem: 


Theorem. Suppose that D and r describe the decrease rate of variances of 


signal x according to Eq. (16). The signal x may be recovered as the solution to 


Eq. (13) with an average value of error 
o cr^iV 


Mt 


, /a^N + MD\ 

‘°8 ( )■ 


(17) 


Equation p7| ) enables us to estimate the number of required samples M of 
signal to achieve a preselected mean recovery error. Intuitively, the af is also 
closely related to the compressibility of signal x, and from Eq. ( |T7| one may 
observe that an average recovery error is inversely proportional to the constant 
value T. The proof of the theorem is given in the Appendix. 
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3. Experimental results 


3.1. Experimental setup 

In this section, we present results illustrating the proposed approach and 
demonstrating that the number of samples (M) required to sense the data can 
be considerably less than the total number of time samples (iV) in the reference 
signal j/°. We investigate the performance of the algorithm using a data set 
of reference signals registered in single module scintillator strip EJ-230 [31] of 
J-PET device [TH] , 

The scheme of the experimental setup is presented in Fig. [l] The 30 cm long 
strip was connected on two sides to the R4998 Hamamatsu |3S| photomultipliers 
denoted as PM1(2). A series of measurements was performed using collimated 
gamma quanta from a ^^Na source placed between the scintillator strip and 
the reference detector. The collimator was located on a dedicated mechanical 
platform allowing it to be shifted along the line parallel to the scintillator strip 
with a submillimeter precision. The ^^Na source was moved from the first to 


PM1 


PM2 


D collimator with 
JJ_I source 


detector 1" 


scope 


Figure 1: Scheme of the experimental setup. 


the second end in steps of 6 mm. At each position, about 5 000 pairs of signals 
from PMl and PM2 were registered in coincidence. These signals were sampled 
using the Serial Data Analyzer (Lecroy SDA6000A) with a probing interval of 
50 ps. To demonstrate the recovery performance only signals from PMl were 
investigated (the procedure with signals from PM2 would be the same). The 
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length of a signal was set to 15 ns, which corresponds to = 300 samples 
(see Fig. [Sjl^for details). 

We wish to make one comment about the data acquisition. The signal cap¬ 
tured by an oscilloscope is length N, where each sample is contaminated with 
white noise with 0 mean and variance. The simulation of measurement y 
is then based on selecting M samples according to the subset fl. However, in 
order to extract the reference, noise-free signal y^, the acquired N samples have 
to be subjected to low pass filtering. In the following procedure we will need 
the signals y and as well. 

Since the absolute registration time has no physical meaning, we synchronize 
the signals in data set Y in such way that the fixed index number 20 corresponds 
to the amplitude of -0.06 V on the rising slope of each signal (see Fig. 

The complete data set Y contains more than 200 000 signal examples and was 
divided into two disjoint subsets: training and testing part, with a ratio 9 to 1, 
respectively. In the training data set only the signals are stored, while in the 
testing one both signals y and are required. 


3.2. Error recovery investigations 

The training data set Y was transformed via PCA into a new space X 
according to the scheme shown in Sec. 2.2. The evaluated matrix A, as well 
as the mean value signal m, were saved and used in the further analysis during 
the signal x recovery from the testing data set. In order to find the theoretical 
value of mean recovery error cr^, introduced in Eq. ( [T7| , one needs to specify 
additionally the following parameters: a,D,T (we will investigate the error 
as a function of the number of samples M). The standard deviation of the noise 
(a) was estimated based on the training data set Y to c.a. 0.015 V, which is 
consistent with the oscilloscope specihcation. The unknown parameters D, t 
were found after the analysis of diagonal elements of the covariance matrix P of 
the training data set X. The smallest value D and the largest value r for which 


the condition from Eq. (16) was met, are equal to 4.2 and 0.33, respectively. 


It should be stressed that, for a given number of samples (M), the expected 
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value of in J-PET scenario would be slightly greater than for the one de¬ 
scribed by Eq. 0- The reason is that in the J-PET scenario the signals are 
probed in the voltage domain and hence in the case when the amplitude of the 
signal is smaller than the threshold level, not all the samples of the signal are 
acquired (see Fig. for example). Therefore, in order to evaluate the theoreti¬ 
cal function of mean recovery error in the J-PET scenario, both, the values of 
the threshold levels as well as the distribution of signal amplitides have to be 
specified first. 

In the first step of the analysis the distribution of signal amplitudes was 
investigated. The experimental cumulative distribution function (cdf), based 
on the signals registered at all the positions along the scintillator strip, is pre¬ 
sented in Fig. The amplitudes of the signals are in the range from -0.3 V to 
-1.0 V. In order to suppress events when gamma quanta were scattered inside 



Figure 2: Experimental cumulative distribution function of signal amplitides. 


the patients body, in the current PET scanners (detecting gamma quanta based 
on photoelectric effect) the energy window, typically in the range from 350 keV 
to 650 keV, is applied [3S] . Such window suppress scattering under angles larger 
than 60 degrees. The J-PET detector is made of plastic scintillators which are 
composed of carbon and hydrogen. Due to the low atomic number of these 
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elements the interaction of gamma quanta with energy of 511 keV is predomi¬ 
nantly due to the Compton effect whereas the interaction via photoelectric effect 
is negligible. In order to suppress scattering in the patient through angles larger 
than 60 degrees, in the J-PET scanner only the signals with energy deposition 
larger than 200 keV will be accepted [S]. Therefore, the signals with amplitude 
smaller than a -0.3 V are filtered out, and a sharp edge of the spectrum for this 
value is seen in Fig. 

In the next step, based on the fully sampled signals stored in testing data 
set Y, we simulate a front-end electronic device that probes the signals at prese¬ 
lected number of voltage levels, both on the rising and falling slope. We carried 
out the experiments for different numbers of voltage levels from 2 to 15. In each 
case, the level of -0.06 V on the rising slope was applied for triggering purposes, 
as was mentioned in Sec. 2.1. The remaining amplitude levels were adjusted 
after a simple optimization process, where the goal was to minimize the exper¬ 
imental mean recovery error tr^. At each step of the optimization process, for 
a fixed number (M) and values of voltage levels, signal recovery was conducted 
in the following way. For each 300 signal samples from testing data set Y, all 
samples at preselected voltage levels were selected to simulate the measurement 
y. Since the amplitude of the signal may be less than certain voltage levels, not 
all samples had to be registered. Therefore, for each processed signal, the num¬ 
ber of acquired samples would be smaller or equal to M. In order to remove the 
mean value from the measurement y, the corresponding values of signal m were 
subtracted from signal samples from the oscilloscope. The measurement matrix 
Aq was formed from the proper rows of matrix A. The signal x was recovered 


using Eq. (13), and finally the signal y was derived as the linear solution of 


Eq. (15). With optimized values of voltage levels, theoretical and experimental 
curves describing the mean recovery error cr^ as a function of the number of 
samples (M) in the J-PET scenario are evaluated and shown in Fig. 

An empirical mean value of is marked with a solid blue line in Fig. 
and is very similar to the expected, theoretical characteristic that takes into 
account the distribution of amplitudes and optimized values of voltage levels 
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Figure 3: Comparison of average recovery errors cr^ as a function of the acquired samples 
(M). Meaning of the curves is described in the text. 


(solid green line). The difference between those two functions is larger for small 
values of M (about 10% of and almost negligible for greater numbers of 
samples. However, both of these functions differ significantly from the theoret¬ 


ical characteristic of cr^, calculated according to Eq.( 17), marked with dashed 
green line in Fig. In the following we will investigate only the case with a 
four-level measurement, which is of most importance since the currently devel¬ 
oped front-end electronic allows one to probe the signals at four fixed-voltage 
levels. It is evident that this comparison of results may be performed in the 
same way for all values of M. 

The optimized values of the four voltage levels are: -0.06, -0.20, -0.35 and - 
0.60 V. Since, the index of the sample taken at the voltage level of -0.06 V at the 
rising slope is common for all signals, the effective number of simulated samples 
at rising and falling edge is equal to M = 7. From Fig. the theoretical value 
of the average recovery error cr^ for M = 7 is c.a. 0.173 (dashed green line). 
However, based on the experimental distribution of amplitudes of the signals, 
presented in Fig. only for about 30% of signals would all samples from four 
thresholds be available (amplitudes larger than -0.60 V). Moreover, for signals 
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with amplitudes in the range from -0.35 V to -0.60 V (about 55% of signals), 
the effective number of samples is equal to 5 and the theoretical value of cr^ 
increases to 0.228 V^. For the rest of the considered signals, with amplitudes in 
the range from -0.30 V to -0.35 V (about 15% of signals), the effective number 
of samples is equal to 3 and the theoretical value of cr^ is 0.346 V^. Finally, the 
expected mean value of cr^ in the J-PET scenario for four voltage levels is equal 
to c.a. 0.227 and is much more comparable with the experimental value 
(equal to c.a. 0.264 V^) than the theoretical value for 7 samples. 

The analysis of the characteristic of tr^ allows us to indicate the proper 
number of samples needed. The function is approximately proportional 


to 1/M but, due to the logaritmic factor (see Eq.( 171), it drops rapidly until 
M reaches the value of about 10. Further increase in the number of samples 
does not provide any significant improvement in the signal recovery. This is 
very important information since the currently developed front-end electronic 
enable one to probe the signals at four fixed-voltage levels, providing eight time 
values for each signal. 

The distribution of the recovery error ||a;° — x\\'^ evaluated using all signals 
from the testing data set for optimized values of four voltage levels is shown 
in Fig. 1^ From the empirical characteristics of ||a;° — i;||| one may see that 
the recovery error is concentrated between 0 and 0.4 with the tail reaching 
the value 1.5 V^. As was shown in Fig. the mean value in the experiment 
is equal to c.a. 0.264 V^. In addition, the standard deviation and the median 
of a probability distribution of a recovery error are equal to c.a. 0.192 and 
0.206 V^, respectively. The three signal recovery examples, with small, medium 
and large recovery error, are shown in Fig. EE and[^ respectively. 

The values of the signal recovery errors in Fig. [^toj^are as follow: 0.082, 
0.266, 0.814 V^. As expected, the worst situation takes a place when the am¬ 
plitude of the signal is slightly below the selected threshold level (see Fig. 
or where it is much larger than the highest sampling voltage. In our sampling 
scheme the highest recovery error occurs for signal amplitudes in the range from 
-0.55 to -0.6 V and from -0.95 to -1 V (where -1 V corresponds to the maxi- 
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Recovery error [V^] 

Figure 4: Distribution of the recovery error evaluated using signals from the testing data set. 



Figure 5: Signal recovery example: the recovery error is about 0.082 V^. 

mum amplitude, see Fig. [^. Unfortunately, there is no possibility to overcome 
these phenomena when only a few samples of the signal are measured. On the 
other hand, it can be seen that the mean value of the error ||a;° — ijlf is on an 
acceptable level. In a typical situation the signal is recovered quite accurately 
(see Fig. [^. 
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Figure 6: Signal recovery example: the recovery error is about 0.266 V^. 



Figure 7: Signal recovery example: the recovery error is about 0.814 V^. 

3.3. Method limitations 

Although, the experimental and theoretical functions describing the recovery 
errors in the J-PET scenario, presented in Fig. are largely consistent, there 
are at least two aspects of the method that need to be investigated: 

1) the assumption about the prior MVN distribution of signals (see 
Eq. which has an impact on the difference bewteen the values of the 


18 























































errors, 

2) the evaluation of the empirical values of a^. as a function of the size of 
training set of signals y®. 

In order to verify the assumption about the normality of signals we have 
used the Kolmogorov-Smirnov test on each of JV principal components in the 
training data set X, evaluated according to Eq. Q. In each dimension, the 
mean value as well as the standard deviation were estimated for all the samples. 
The significance level used in this study was 0.05. The hypothesis, regarding 
the normal distribution form, was rejected only for the first principal component 
that holds about 40% of the signal energy. However, in that case the calculated 
value from the statistical test was not significantly higher than the critical value. 
From this analysis one infers that the signals stored in matrix X are not exactly 
normally distributed. This fact may contribute to the difference between the 
theoretical and empirical values of cr^. 

It should be stressed that all the informations needed to recover the signal x 
(matrix A, P and vector /r - see Sec. 2.2 for details), and therefore the empirical 
values of were evaluated based on large set of about 200 000 training signals 
y°. In the following we will analyze the influence of the size of the training set 
of signals on the value of the signal recovery error. We conduct the experiment 
for a wide range of number of signals y° in the training set from 50 to 200 000. 
In each case we investigate only the four-level measurement y [M = 7). The 
results of the analysis of the empirical values of cr^ as a function of the size of 
the training set of signals y° are shown in Fig. ([^. 

The parameters of the model that were used for the signal recovery in the 
study described in Sec. 3.2 were based on 200 000 training signals, which corre¬ 
sponds to the value of the empirical recovery error of 0.264 V^. From Fig. (|^ 
one may observe that reducing the number of training signals down to about 
10 000 does not influence the quality of the recovery of signal i; the cr^ error 
is almost constant in that range. However, for smaller number of traning sig¬ 
nals, the error increases rapidly and the recovery of the signal x becomes 
increasingly less accurate. 
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Figure 8: Influence of the size of the training set of signals on the average recovery error 


3.4- Spatial resolution of the hit-position reconstruction 

In this section, we will incorporate the method for hit-position reconstruc¬ 
tion, described in Ref. [15], in order to evaluate a position resolution of the 
J-PET scanner with fully recovered signals. We will compare the spatial reso¬ 
lutions obtained from the original raw-signal (300 samples) to those from the 
compressed signal (e.g. 8 samples). We have carried out experiments with num¬ 
bers of voltage levels from 2 to 15, which corresponds to the number of samples 
M from 3 to 29. 

For a single event of gamma quantum interaction along the scintillator strip, 
a pair of signals at two photomultipliers is measured in a voltage domain. Next, 
the signals are recovered according to the description in Sec. 2, and finally, an 
event is represented by a 600-dimensional vector. For a fixed number of voltage 
levels a two-step procedure of the position reconstruction was performed. First, 
the scintillator’s volume was discretized and for each bin a high statistics set of 
reference 600-dimensional vectors was created. The objective of the second part 
of the procedure is to classify the new event to one of the given sets and hence 
determine the hit position. For more details about conducting the experiment 
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of hit position reconstruction, the reader is referred to Ref. m- 

We have conducted the test on the same data set and under the same con¬ 
ditions as described in Ref. HD, where the spatial resolution was reported to be 
equal to 1.05 cm (a). The spatial resolutions derived from the recovered signals 
as a function of the number of samples M included in the recovery process are 
shown in Fig. 



Figure 9: The spatial resolution as a function of the acquired samples (M). 


In Fig. only the region for small M, from 3 to 29, is shown, but it has to 
be stressed that the spatial resolution derived from the original raw-signal (300 
samples) is equal to 0.933 cm {a) and is almost the same as for M=29. For the 
most interesting case, with four voltage levels (M=7), the spatial resolution is 
slightly worse than for the fully sampled signal and is equal to 0.943 cm (cr). On 
the other hand, even in that case the spatial resolution is about 0.1 cm better in 
comparison to the one evaluated based on signals in the voltage domain alone. 

4. Conclusions 

In this paper a novel scheme of recovery of signals generated in plastic scintil¬ 
lator detectors in the J-PET scanner was introduced. The idea of signal recovery 
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is based on the Tikhonov regularization theory, that uses the training data set 
of signals. In these studies we assumed that training signals come from a MVN 
distribution. The compact representation of these signals was provided by the 
PCA decomposition. 

One of the most important aspect of our work considers a statistical analysis 
of an error level of recovered signals. In this work a dependence of the signal 
recovery error on the number of samples taken in the voltage domain was de¬ 
termined. It has been proven that an average recovery error is approximately 
inversely proportional to the number of samples and inversely proportional to 
the decrease rate of variances in the covariance matrix. In the experimental 
section, the method was tested using signals registered by means of the single 
detection module of the J-PET detector. It was shown that the PCA basis 
offers high level of information compression and an accurate recovery may be 
achieved with just 8 samples for each signal waveform. It is worth noting that 
the developed recovery scheme is general and may be incorporated in any other 
investigation where a prior knowledge about the signals of interest may be uti¬ 
lized. 

In the experimental section we have demonstrated that using the recovered 
signals improves the hit-position reconstruction. In order to evaluate a position 
resolution of the J-PET scanner with fully recovered signals, we have incorpo¬ 
rated the method for hit-position reconstruction, described in Ref. m- In the 
cited work, the spatial resolution evaluated on the same data set and under 
the same conditions, based on 8 samples in voltage domain, without a recovery 
of the waveform of the signal, was reported to be equal to about 1.05 cm {a). 
Our experiment shows that the application of an information from four voltage 
levels to the recovery of the signal waveform can improve the spatial resolution 
to about 0.94 cm (cr). Moreover, the obtained result is only slightly worse than 
the one evaluated based on all 300 samples of the signals waveform. The spatial 
resolution calculated under these conditions is equal to about 0.93 cm {a). It is 
very important information since, limiting the number of threshold levels in the 
electronic devices to four, leads to a reduction in the cost of the PET scanner. 
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Future work will address a development of the more advanced method to 
define the hit-position and event time for annihilation quanta in the J-PET de¬ 
tector based on the recovered information. We believe that, with fully recovered 
signals, there is still scope for improvement in the time and position resolution 
of the J-PET scanner. 


Appendices 


In order to prove the theorem we assume, for the sake of simplicity, that the 
matrix A has normally distributed elements with zero means and 1 /IV variances. 
These values of the parameters of normal distribution ensure that the matrix A 


is orthonormal. Hence, based on Eq. (141, the matrix S is given by: 

-1 


5'= P 


3-1 




The cr^ is equal to the trace of the matrix S and hence: 


N 


(18) 


= E 


a^NPk,k 
am + MPk,k 


2^2 


N 


M 




1 


(19) 


k=l 


am + MPkk 


The sum in the last term in Eq. (19) may be approximated by a definite integral. 
In the following we will use for the calculations a basic, rectangle rule, and: 

N . pN+h 


E 

fc=i 


1 


1 


a^N + MPk,k Ji-h am + MP{k) 


dk = l 


where h = 1/2. At the very beginning we assumed that the function P{k) has 
the form: P{k) = D ■ e~'^^ (see Eq. @)- We will perform the integration using 
the substitution t = e~'^^. Without any significant loss of precision, we change 
the integration limits from [1 — li, N + h] to [0, N], The calculations of the 
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integral I will be as follow: 


J^-^N {a"^N + MDt) 
/*! 1 />! 
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g-TW a’^Nrt 
1 


dt — 


MD 


a'^Nr 


-tn a^Nr {a'^N + MDt) 
(log(t)|g_^N - log(cr^iV + MDt)\],-rN) 


dt 
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t'^Nt 


Nt + log 


cf^N 


a'^N + MD 


and thus 
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